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We report on a measurement of the cosmic ray energy spectrum with the IceTop air shower array, the 
surface component of the IceCube Neutrino Observatory at the South Pole. The data used in this analysis 
were taken between June and October, 2007, with 26 surface stations operational at that time, corre- 
sponding to about one third of the final array. The fiducial area used in this analysis was 0.122 km * 1 2 3 4 . 
The analysis investigated the energy spectrum from 1 to 100 PeV measured for three different zenith 
angle ranges between 0° and 46°. Because of the isotropy of cosmic rays in this energy range the spectra 
from all zenith angle intervals have to agree. The cosmic-ray energy spectrum was determined under dif- 
ferent assumptions on the primary mass composition. Good agreement of spectra in the three zenith 
angle ranges was found for the assumption of pure proton and a simple two-component model. For 
zenith angles 9 < 30°, where the mass dependence is smallest, the knee in the cosmic ray energy spectrum 
was observed at about 4 PeV, with a spectral index above the knee of about -3.1. Moreover, an indication 
of a flattening of the spectrum above 22 PeV was observed. 

© 2013 Elsevier B.V. All rights reserved. 


1. Introduction 

100 years after the discovery of cosmic rays, their sources and 
acceleration mechanisms still remain mostly unknown. The energy 
spectrum of cosmic rays as measured by various experiments fol- 
lows a relatively smooth power law with spectral index y « -2.7 
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up to about 4 PeV, where it steepens to y « -3.1 [1]. While this 
feature in the spectrum called “knee” is well established, its origin 
remains controversial [2]. Most models to explain the knee involve 
a change in chemical composition of cosmic rays in the energy 
region above the knee. Such a change has been observed by various 
experiments [3] but systematic uncertainties are too large to 
discriminate individual descriptions. Features in the all-particle 
cosmic ray energy spectrum and their chemical composition bear 
important information on the acceleration and propagation of 
cosmic rays. The measurement of the cosmic ray energy spectrum 
and composition is the main goal of the IceTop air shower array. 

IceTop is the surface component of the IceCube Neutrino Obser- 
vatory at the geographic South Pole [4]. Installation of IceCube and 
IceTop was completed at the end of 2010, with 86 IceCube strings 
and 81 IceTop stations deployed covering an area of about 1 km 2 
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and a volume of about 1 km 3 . IceTop was designed to measure the 
energy spectrum and the primary mass composition of cosmic ray 
air showers in the energy range between 5 ■ 10 14 eV and 10 18 eV. 

The average atmospheric depth at the South Pole is about 680 g / 
cm 2 . IceTop is therefore located close to the shower maximum for 
showers in the PeV range (for vertical protons about 550 g/cm 2 at 
1 PeV to 720 g/cm 2 at 1 EeV). This has the advantage that local 
shower density fluctuations are smaller than at later stages of 
shower development. 

In this paper, we present the first analysis of IceTop data on 
high-energy cosmic rays and a measurement of the cosmic ray en- 
ergy spectrum. This analysis is based on air shower data taken with 
the IceTop surface stations. The data were taken between June and 
October 2007 with 26 IceTop stations operating, which comprise 
about 1 /3 of the complete detector. 

Section 2 of this paper gives an overview over the IceTop array, 
and the processing and calibration of tank signals, which are the 
basis for reconstructing air showers. Section 3 describes the data- 
set and run selection criteria. Section 4 introduces event recon- 
struction, and in Section 5, simulation of air showers and of the 
IceTop tank response are presented. In Section 6 the final event 
selection and detector performance are discussed. Section 7 de- 
scribes the determination of the primary energy, whereas system- 
atic uncertainties are discussed in Section 8. In Section 9 the results 
are presented and discussed. 

2. The detector 

2.1. IceTop 

The IceTop air shower array is the surface component of Ice- 
Cube, covering an area of about 1 km 2 with 81 detector stations 
above the 86 IceCube strings [4]. The stations are mostly located 
next to IceCube strings with a average spacing of 125 m, except 
for three stations placed as an infill with a smaller spacing in the 
central part of the detector, in order to lower the energy threshold 
of the detector to about 300 TeV. By 2007, 22 IceCube strings and 
26 IceTop stations had been deployed. These stations are 


highlighted in Fig. 1, which shows the layout of the IceTop air 
shower array in its final configuration. 

Each station consists of two ice-filled tanks separated from each 
other by 10 m. The two tanks of each station are embedded in 
snow with their tops aligned with the surface in order to minimize 
the accumulation of drifting snow (see Section 2.5) and to protect 
the ice from temperature variations. 

The tanks are cylindrical with an inner diameter of 1 .82 m, and 
are filled with transparent ice to a depth of 90 cm (see Fig. 2). The 
inner tank walls are covered with a diffusely reflective coating. The 
first four stations deployed in 2005 have a liner with a higher 
reflectivity. This difference affects amplitude and pulse width of 
detected tank signals, since the higher reflectivity reduces Cheren- 
kov photon absorption, leading to longer pulses. The ice is covered 
with perlite as thermal and optical insulation. The perlite also of- 
fers proper diffuse reflectivity at the ice surface. 

Each tank is equipped with two ‘Digital Optical Modules’ 
(DOMs) [5] to record Cherenkov light generated by charged parti- 
cles passing through the tank. The DOMs are identical to those 
used in other IceCube components and consist of a 10" photomul- 
tiplier tube (PMT) [6], plus electronic circuitry for signal digitiza- 
tion, readout, triggering, calibration, data transfer and various 
control functions. The two DOMs in each tank were operated at dif- 
ferent PMT gains, 5 ■ 10 6 (high-gain DOM) and 5 • 10 5 (low-gain 
DOM), to enhance the dynamic range. This resulted in a linear dy- 
namic range from 1 to more than 10 5 photoelectrons (PE). During 
the data taking period used in this analysis all 104 DOMs in the 26 
IceTop stations were fully operational. 

2.2. Trigger and data acquisition 

A DOM records PMT signals autonomously. A signal is recorded 
if it surpasses a certain discriminator threshold, which in the case 
of IceTop was set to 22 mV for the high-gain DOMs (corresponding 
to about 20 pe) and 12 mV for the low-gain DOMs (corresponding 
to about 180 pe). The exact charge threshold depends on the shape 
of the IceTop multi-photoelectron pulses, which is determined by 
the arrival times of photoelectrons. After triggering, the delayed 
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Fig. 1. Layout of the IceTop air shower array. Colors indicate the year of deployment and the 26 stations installed in 2007 are highlighted. 
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Fig. 2. Cross section of a tank showing the tank geometry with insulation and position of the DOMs. The center of the ice surface between the two DOMs is used as tank 
position by reconstruction algorithms. 



PMT pulse is sampled by ‘Analog Transient Waveform Digitizers’ 
(ATWDs) with three different gain channels (nominal gains are 
0.25, 2, and 16) in 128 bins with a width of 3.3 ns, corresponding 
to a total sampling time of about 422 ns. The analog samples are 
then digitized to 10 bit accuracy. 

Up to this point, signal recording happens independently in 
each DOM. To reduce the high trigger rates in high-gain DOMs 
(~2 kHz), which are mostly from low-energy showers, a hardware 
‘local coincidence’ with a coincidence time window of ±1 ps be- 
tween the high-gain DOMs in the two tanks of a station is required 
to initiate the readout and transmission of DOM data to the count- 
ing house (IceCube Lab). 

In the counting house an event is built if the IceTop software 
trigger condition is satisfied requiring six or more DOMs to report 
a (local coincident) signal within a time window of 5 ps. The event 
will contain waveforms from all IceTop and IceCube DOMs, which 
fulfilled the local coincidence condition between 10 ps before the 
first until 10 ps after the last of the six DOMs that triggered the 
event building. The requirement of 6 DOMs means that at least 
two stations had to trigger. For the 26 station configuration, the to- 
tal IceTop trigger rate was about 14 Hz. 

2.3. Charge extraction and calibration 

Fig. 3 shows a typical waveform measured in IceTop. While 
waveforms are recorded in three ATWD channels, this analysis 
used only the highest gain unsaturated (less than 1022 ADC 
counts) channel. In this analysis only the integrated charge (i.e. 
the integral over the whole 422 ns waveform) and the signal time 
were used. Before a waveform was integrated, its baseline was sub- 
tracted by determining the average value in bins 83 to 123 high- 
lighted in the figure. The undershoot is caused by droop 
introduced by the transformer used to couple the photomultiplier 
tube to the DOM’s front-end electronics. The signal time (‘leading 
edge time’) was defined by extrapolating the steepest rise of the 
waveform before the maximum down to the baseline. The absolute 
time scale of a DOM is calibrated with respect to all other DOMs to 
an accuracy of about 3 ns [7]. 

The charge produced by a single photoelectron, the amplifier 
gains and the digitizers are calibrated in a procedure common to 
all IceCube DOMs [7]. However, the signal response to a particle 


of a given type and energy traversing the tank, expressed in photo- 
electrons, differs from tank to tank, due to differences in ice quality 
and reflectivity of the tank walls. Therefore, the signal of each tank 
is converted to a common unit called ‘Vertical Equivalent Muon’ 
(VEM). Calibration was done by recording charge spectra of DOMs 
in dedicated calibration runs with all DOMs operated at a gain of 
5 10 6 and without requiring local coincidence (for an example 
see Fig. 3, right). These charge spectra show a clear peak due to pe- 
netrating muons above a background of electrons and photons. The 
spectra are fitted by the sum of a function describing the muon 
peak and an exponentially falling background term. Measurements 
with a portable scintillator telescope mounted on top of tanks, 
restricting muons to nearly vertical angles of incidence, indicated 
that the peak for vertical muons lies about 5% lower than for the 
full angular range. Simulation studies confirmed that restricting 
the angles of incidence of muons shifts the peak position by about 
5% [8], The scaled peak is referred to as ‘VEM peak’. For a given 
DOM the VEM unit can be expressed in terms of number of photo- 
electrons. These values average 120 and 200 photoelectrons for the 
low and high reflectivity tanks (see above), respectively. 

For the 5-month run, 15 calibration runs were used. Between 
two consecutive calibration runs, the charge calibration was as- 
sumed to be stable (see also the discussion in Section 8.4). 

2.4. Atmospheric conditions 

Variations of the atmosphere influence the development of air 
showers and thus the signals measured in IceTop. Since IceTop is 
below the shower maximum for all energies of interest in this anal- 
ysis and for all primary masses, an increase of the atmospheric 
overburden leads to an attenuation of shower sizes. Vertical atmo- 
spheric overburden is related to ground pressure p as X v = p/g, 
where g= 9.87 m/s 2 is the gravitational acceleration at the South 
Pole. While there is some annual variation of the ground pressure, 
it mostly varies on shorter time scales on the order of days. 

Besides ground pressure, the altitude profile of the atmosphere, 
dX„(h)/dh, also influences the development of air showers. This 
altitude profile has a pronounced annual cycle because the cold 
atmosphere during the winter months is much denser than the 
warmer atmosphere of the summer months. The data used in this 
analysis were mostly taken during the winter months. 
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Fig. 3. Left: A typical IceTop waveform. The blue horizontal line marks the baseline and the near vertical green line indicates the extrapolation of the leading edge yielding the 
signal time marked by the red circle. The baseline is below 0 (dashed line) due to droop. Right: A typical charge spectrum recorded for the VEM calibration. The spectra are 
fitted with an empirical formula to determine the peak position. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of 
this article.) 


In the simulations used to interpret the air shower data a model 
of the South Pole atmosphere is used, which should represent the 
average atmosphere during the data taking period. Nevertheless, 
variations of the atmosphere around the average lead to an addi- 
tional uncertainty on the measured energy spectrum. These sys- 
tematic uncertainties will be discussed in Section 8.2. 


2.5. Snow 

During installation, IceTop tanks are embedded in snow up to 
the upper surface of the tanks. Depending on location, surrounding 
surface and structures, each tank is covered by accumulated layers 
of snow of varying thickness. Each year the amount of snow on the 
IceTop tanks grows on average by 20 cm. 

As shown in Fig. 4, the snow height for the analyzed data varied 
mostly between 0 and 30 cm, except for four stations close to a 
building, which are covered by 60 to 90 cm of snow. The average 
snow height was 20.5 cm in January, 2007. 

The snow has an average density of 0.38 g/cm 3 , depending on 
snow height and location. The snow on top of and around the tanks 
influences the response to air shower particles penetrating the 
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Fig. 4. Snow heights on top of IceTop tanks measured in January 2007. All 20 newly 
deployed tanks had no snow on top, the average snow height was 20.5 cm. The 
dashed histogram is the snow height distribution on top of the same tanks 
measured one year later. 


tanks and needs to be taken into account in simulations and for 
the determination of the shower energy. 


3. Data set and data selection 

Event filtering and data transmission. The data used in this anal- 
ysis were taken between June 1st and October 31st, 2007. The anal- 
ysis was performed using a data sample which was transferred 
with limited bandwith via satellite to the IceCube data center at 
UW Madison. Due to these bandwidth constraints, events with less 
than 16 participating DOMs were prescaled by a factor of 5. Events 
with 16 or more DOMs were transmitted at a rate of 0.9 Hz and the 
small events at a rate of 2.5 Hz. Due to the high event rate, the pre- 
scale does not cause any statistical problems in the low-energy 
region. 

Run selection. In order to ensure detector stability and data qual- 
ity, the following criteria were applied to runs which were used in 
this analysis: 

• The run was longer than 30 min. A normal detector run lasted 
8 h, and nearly all runs that were aborted after a short time 
encountered some sort of problem. 

• All DOMs were running stably. 

• After correction for atmospheric pressure variations the trigger 
and filter rates were stable and within ±5% agreement with the 
previous good run. Pressure correction was done by fitting the 
relation between ground pressure p and rate R with an expo- 
nential function, R(p) ~ exp (-bp), yielding a barometric coeffi- 
cient b = 0.0077/mbar [9]. Then, the rates were corrected to 
the average South Pole ground pressure of 680 mbar: 


^corrected = Rexp(b(p - 680 mbar)). (1) 

These cuts reduced the livetime by about 10%. 

Event cleaning. Before starting the reconstruction, events were 
cleaned based on a few simple timing criteria. When both DOMs 
of a tank triggered, the tank signal was rejected if the time differ- 
ence between the two signals was greater than 40 ns. The analysis 
used only one signal per tank. For each high-gain DOM a saturation 
threshold was determined from a comparison of signals that trig- 
gered both DOMs in a tank. Signals with less charge were taken 
from the high-gain DOM. If the charge exceeded the saturation 
threshold, the charge measured by the low-gain DOM was 
used and the time was determined from the high-gain signal. 
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Furthermore, a tank signal was also rejected if only the low-gain 
DOM triggered and the high-gain DOM was missing. 

Then, a maximum time difference of 

\t A -t B \ < |Xyl ~ Xg| + 200 ns (2) 

between signals in tanks A and B of the same station was required. 
Here, t A and t B are the signal times in the two tanks and x A and x B 
are the tank locations. The tolerance of 200 ns was introduced in or- 
der to account for shower fluctuations. Finally, stations were 
grouped in clusters, such that any pair of stations i and j in the clus- 
ter fulfilled the condition 

\U - tj| < |Xi ~ Xjl + 200 ns. (3) 

The station position x t is the center of the line connecting its two 
tanks, and t, is the average time of the tank signals. In each event, 
only the largest cluster of stations was kept. 

Only about 10% of events were affected by this event cleaning, 
and about 2.5% of events dropped below the threshold of 5 stations 
required for reconstruction. On average 2.3 tanks were removed. 

Charge-based retriggering. In order to reduce uncertainties due 
to the description of the detector threshold in the simulation, all 
events were retriggered to a common threshold based on total reg- 
istered charge. All pulses with a charge below S t h r = 0.3 VEM were 
removed, and afterwards the local coincidence conditions (see Sec- 
tion 2.2) were re-evaluated discarding all pulses that no longer ful- 
filled this condition. This procedure was applied to both 
experimental and simulated data. 

Event selection. For further processing, a total of N tot = 
8,895,205 events were selected where at least five stations had 
triggered. These five stations do not have to be neighbors. Events 
which fulfilled this condition, but had less than 16 DOMs read 
out (before event cleaning), were reweighted in the analysis with 
the prescale factor of 5 (see above). 

The effective livetime was calculated by fitting the distribution 
of time differences between events, At, with an exponential 
function, 

N( At) = N 0 exp(-At/t). (4) 

This was done individually for each data taking run. The selected 
runs have a total effective livetime of T = JjmnsiW ■ t,j = 
(3274.0 ± 1.9) h, which corresponds to 89.4% of the selected 
153 days period. The uncertainty on the livetime was included in 
the statistical error. 

4. Air shower reconstruction 

The energy of the primary particle cannot be measured directly, 
but has to be determined from the properties of the observed air 
shower. From the measured charges and times of the tank signals 
the shower core position, the shower direction, and the shower 
size are reconstructed. The latter is a measure of primary energy 
and was defined here as the signal S125 measured at a distance 
R = 125 m from the shower axis. Details of the reconstruction pro- 
cedure are described in [4,10]. In the following we summarize the 
essential steps. 

Air shower reconstruction required at least 5 triggered stations. 
The logarithm of the signal, logS(R), of a tank at distance R from the 
shower axis was fitted by a lateral distribution function which is a 
second order polynomial in logR [11]: 

logS(R) = logS, 25 - ft log (1 J - Klog 2 (125 m ) ■ 0 ) 

The free parameters of the function, in addition to the shower size, 
S125, are the slope j) and the curvature K. The latter parameter was 


a 




Fig. 5. Left: Example of an IceTop lateral fit. The shower triggered 25 stations and 
the reconstructed shower size is S 12 s = (65.1 ± 2.8) VEM. Right: Time residuals with 
respect to a plane perpendicular to the shower direction fitted by the sum of a 
parabola and a Gaussian function. "Upstream" and “downstream" refer to tanks 
being hit before and after the shower core reaches the ground. 


fixed to k = 0.303, the average value found in simulation studies. 
Implicitly the function (5) depends on four more parameters 
since R depends on the shower core position (x c ,y c ) and direction 
(0, <j>). Since the function (5) behaves unphysically at small dis- 
tances to the shower axis (R < 1 m) all signals within 1 1 m of 
the core are excluded from the fit. Fig. 5 left shows an example 
of the lateral distribution function fit of a shower with 25 trig- 
gered stations. 

The arrival times of the signals map out the shower front which 
was described by the sum of a parabola and a Gaussian function, 
both symmetric around the shower axis (Fig. 5 right). 

The lateral distribution function and the function describing the 
shower front curvature the charge and time of air shower signals 
were fitted to the measured signal charges and times using the 
maximum likelihood method. In addition to terms for the signal 
charges and times, the likelihood function also takes into account 
stations that did not trigger. 

5. Simulation of air showers and the IceTop detector 

The relation between the measured signals and the energy of 
the primary particle, as well as detection efficiency and energy res- 
olution were obtained from CORS1KA [12] air shower simulations 
and simulations of the IceTop detector. 
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5. 1. Air shower simulation 

We simulate the development of air showers in the atmosphere 
using the simulation code CORSIKA [12]. Inside CORSIICA, the ha- 
dronic component of the air showers was simulated using the 
models S1BYLL2.1 [13,14] and FLUKA 2008.3 [15,16] for the high 
and low energy interactions, respectively. The electromagnetic 
component was simulated using the EGS4 code [17] and no ‘thin- 
ning’ (reduction of the number of traced particles) was applied. To 
study systematic effects of the hadronic interaction model, small 
samples of showers were simulated using the QGSJET-11 [18,19] 
and EPOS 1.99 [20] high energy interaction models. Two different 
parameterizations of the South Pole atmosphere from two days 
in 1997 based on the MSIS-90-E model [21] were used: July 1st 
and October 1st (CORSIKA atmospheres 12 and 13). The July atmo- 
sphere has a total overburden of 692.9 g/cm 2 , while the October 
atmosphere has an overburden of 704.4 g/cm 2 . The July atmo- 
sphere was used in the data analysis, because its total overburden 
is close to the average measured overburden of 695.5 g/cm 2 and its 
profile corresponds to that of a South Pole winter atmosphere. The 
October atmosphere model was used to study systematic uncer- 
tainties due to the atmospheric profile used in the simulation. 

5.2. Detector simulation 

The output of the CORSIKA program, i.e. the shower particle 
types, positions and momenta at the observation level of 2835 m, 
were injected into the IceTop detector simulation. The simulation 
determines the amount of light produced by the shower particles 
in the tanks followed by the simulation of the PMT, the DOM elec- 
tronics and the trigger chain. 

The Cherenkov emission inside the tanks is simulated using 
Geant4 [22,23]. All structures of the tank, the surrounding snow, 
including individual snow heights on top of each tank, as well as 
the air above the snow are modeled realistically [4]. The snow 
heights used in the simulation corresponded to those measured 
in January 2007 (see Fig. 4). In order to save computing time, Cher- 
enkov photons are not tracked; only the number of photons emit- 
ted in the wavelength interval 300 nm to 650 nm is recorded. 
Using Geant4 simulations, that include Cherenkov photon tracking 
until photons reach the PMT, it was shown that the number of de- 
tected photons scales linearly with the number of emitted photons, 
independent of incident particle type and energy. The propagation 
of Cherenkov photons is modeled by distributing the arrival times 
according to an exponential distribution, which is tuned such that 
simulated waveform decay times match those observed in experi- 
mental data (26.5 ns for all tanks commissioned after 2005 and 
42.0 ns for tanks commissioned in 2005). 

The number of photoelectrons corresponding to 1 VEM was ta- 
ken from the VEM calibration of the real tanks and used as an input 
for the simulation. The simulated tanks were then calibrated by 
generating muon spectra as in experimental data using air shower 
simulations with primary energies between 3 GeV and 30 TeV and 
zenith angles up to 65°. Thus, the ratio between the number of 
emitted Cherenkov photons and observed photoelectrons was 
determined by the VEM calibration of simulated tanks. 

In the next step the generated photoelectrons are injected into a 
detailed simulation of the PMT followed by the analog and digital 
electronics of the DOM. To simulate the photomultipliers, Gaussian 
single photoelectron waveforms with a random charge according 
to the average single photoelectron spectrum are superimposed 
[6]. Afterwards, a saturation function is applied to the resulting 
waveforms. In the DOM simulation, the pulse shaping due to the 
analog front end electronics is applied to the output of the PMT 
simulation. This includes the individual shaping of the signal paths 
to the ATWD and the discriminators, as well as the simulation of 


the droop effect induced by the toroid that couples the high voltage 
circuits of the PMT to the readout electronics. Then, the discrimina- 
tors are simulated and the local coincidence conditions are evalu- 
ated. Finally, the waveform digitization and the array trigger are 
simulated. 

Simulated data are of the same format as the experimental data 
and were reconstructed in the same way, as described in the pre- 
vious section. 


5.3. Simulation datasets 


In this analysis we describe the cosmic ray composition just 
with the two extreme nuclei proton and iron. The justification 
comes from the fact that the final result is not very sensitive to de- 
tails of the composition assumption but mostly to the mean loga- 
rithmic mass. 

In total 2 • 10 5 showers of proton and iron primaries in the en- 
ergy range between 100 TeV and 100 PeV were generated in 30 
logarithmic energy bins according to an E 3 distribution. For the 
analysis, the events are reweighted to an E~ 3 flux, which is closer 
to the results of previous experiments and thus reduces systematic 
biases (see also Section 8.10). In addition to pure proton and iron 
simulations we also combined the datasets using a parametriza- 
tion of Glasstetter’s two-component model [24] (see Fig. 6). We 
transformed the proton flux to the form 


d / / e y> +1 

dE ° VlPevJ 



( 6 ) 


as suggested in [25], with f 0 = 3.89 ■ 1CT 6 m _2 s _1 sr _1 , 
= -2.67, y 2 = -3.39,£|< n ee = 4.1 PeV, and e = 2.1. The iron flux 
was used as specified in Ref. [24]: 


^ = 1.95 ■ 10~ 6 itt 2 s 
dt 


1 PeV 


( 7 ) 


The total flux was then normalized to the same E 3 spectrum as in 
case of the single component Monte Carlo. 

Since shower generation is CPU intensive the same showers 
were sampled several times inside a circle with a radius of 



log(E 0 /PeV) 


Fig. 6. Relative abundance of proton and iron in our parametrization of Glasstetter’s 
two-component model as a function of primary energy. Above about 10 PeV the 
spectrum is dominated by iron. 
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1200 m around the center of the 26 station IceTop array. The num- 
ber of samples was chosen for different energy bins such that every 
shower would remain on average only once in the final sample 
after applying the cuts described in the next section. This ensures 
a good balance between an effective use of the generated showers 
and the artificial fluctuations introduced by oversampling. 

6. Event selection and reconstruction performance 

Quality cuts. Based on the reconstruction results the following 
quality criteria were required for each event entering the final 
event sample, for both simulated and experimental data: 

• Containment cut: The reconstructed core and the first-guess 
core position had to be at least 50 m inside the boundary of 
the array. The array boundary is defined by the polygon with 
vertices at the centers of stations at the periphery of the array 
and edges connecting these stations. This cut defines a fiducial 
area of A cut = 0.122 km 2 . Furthermore, it was required that the 
station containing the largest signal is not on the border of 
the array. 

• Only events with zenith angles 0 < 46° were considered. 

• The reco nstructio n uncertainty on the core position had to fulfill 
(T C ore = \J< 7* + aj < 20 m, which is the fit uncertainty of the 
position parameters calculated from the width of the likelihood 
function around the minimum. 

• The slope parameter /i had to be in the range 2.0 ^ < 4.5 

because most events with fi values outside this range were 
badly reconstructed and because j) was limited in the fit. The 
eliminated events predominantly had low primary energies, 
E 0 < 1 PeV. 

In the experimental dataset 3,096,334 events passed the quality 
cuts. Passing rates for the individual cuts are shown in Table 1 for 
events with S125 > 1VEM. Differences between data and two-com- 
ponent Monte Carlo are discussed later in Section 8.7. 

Reconstruction performance. Core position and angular resolu- 
tion, shown in Fig. 7, are key criteria for the performance of air 
shower reconstruction. The lcr core resolution is defined as the 
68% quantile of the cumulative distribution of the distances 
between true and reconstructed shower cores; correspondingly 
the angular resolution is defined as the angle between true and 
reconstructed shower direction. The numbers shown are for show- 
ers with zenith angle 0 r' : 30°, obtained from the two-component 
Monte Carlo after applying the quality cuts listed in the previous 
paragraph. At the highest energies, a core resolution of 7 m and 
an angular resolution of 0.4° were achieved. In the most inclined 
zenith angle range considered in this analysis, 40° ^ 9 < 46°, the 
core resolution is between 10 and 30% worse and the angular 
resolution between 10 and 25% worse, depending on energy and 
primary mass. 


Table 1 

Passing rates of the quality cuts described in the text for events with S , 2 5 > 1 VEM. 
Statistical errors on experimental data are negligible. 


Cut 

Experimental data 

Monte Carlo 



Passing Cumulative 
rate (%) (%) 

Passing rate 

(%) 

Cumulative 

(%) 

^station > 5 and 

100 


100 


S 125 > 1 VEM 

Largest signal contained 

42.5 

42.5 

(39.4 ±0.5) 

(39.4 ±0.5) 

First guess core 

95.8 

40.7 

(95.4 ±0.4) 

(37.6 ±0.5) 

contained 

Core contained 

78.9 

32.1 

(81.1 ±0.5) 

(30.5 ±0.6) 

Zenith 0 < 46° 

96.3 

30.9 

(96.4 ±0.5) 

(29.4 ±0.6) 

CT C ore 20 m 

99.7 

30.8 

100 

(29.4 ±0.6) 

2.0 < j6 < 4.5 

98.1 

30.2 

(99.7 ±0.1) 

(29.3 ±0.6) 




Fig. 7. Core position and angular resolution for showers with 0 c: 30°, obtained 
from the two-component Monte Carlo. At high energies, the distance between true 
and reconstructed core position of 68% of showers is 7 m or less. The angle between 
true and reconstructed directions of 68% of showers at 1 PeV is smaller than 0.8° 
and this value decreases to 0.4° at 100 PeV. 


In the most inclined zenith angle range considered in this anal- 
ysis, 40° ^ 6 < 46°, a core resolution of 10 m and an angular reso- 
lution of 0.5° was achieved. 

The intrinsic spread of the shower size distribution due to the 
detector was studied using such simulated proton air showers, 
which remained more than once in the final sample after resam- 
pling, detector simulation, and reconstruction. While the value of 
S 125 cannot be predicted or calculated analytically for a given air 
shower, the distribution of reconstructed shower sizes for a given 
air shower (at different locations inside the detector), is a measure 
of the intrinsic resolution of the apparatus and the measurement. 
As shown in Fig. 8, the RMS of log 10 (Si 25 /V£M) improves from 
about 0.08 at 1 PeV primary energy to better than 0.04 above pri- 
mary energies of 1 0 PeV, almost independent of zenith angle. A 
comparison of this intrinsic resolution with the spread of recon- 
structed energies discussed in Section 7.2 shows that for the most 
vertical zenith angle range {0 < 30°) the intrinsic resolution is 
dominant, whereas the energy resolution of more inclined showers 
is dominated by fluctuations of the shower development. Recon- 
structed shower core distribution. The distribution of recon- 
structed shower core locations is shown in Fig. 9. While this 
distribution is not entirely flat, it is reasonably well reproduced 
by simulation, and the structures are understood as described be- 
low. Therefore, we do not expect a significant systematic error 
introduced by the containment criteria. 
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Fig. 8. RMS spread of the S , 2 5 distribution due to the apparatus and the shower 
reconstruction plotted against the primary energy. 



X / m 


Fig. 9. Distribution of reconstructed core positions in IceTop for showers with 
Si 25 3 s 2 VEM. The black line represents the containment region. The dots are the 
IceTop tanks, red open circles being those that are considered on the border of the 
array for the containment cut. (For interpretation of the references to colour in this 
figure legend, the reader is referred to the web version of this article.) 

The l ight half of the array (x > 200 m) has been deployed in 
2005/06 and the left half has been deployed in 2007. Consequently, 
the right half of the array is covered by a thicker layer of snow. This 
increases the energy threshold in this part of the array, leading to a 
lower rate for cores with x > 200 m. 

There are regions along the border of the containment region 
(black line), where the closest station is one that is not considered 
contained (red dots in the figure). Due to the requirement that the 
station with the largest signal should be contained, shower cores 
reconstructed in this area are less likely to pass the containment 
cut. It should be noted that these stations still contribute to the lat- 
eral fits of all events. 

Finally, there are some structures in the vicinity of stations, 
which are an artifact of the removal of tanks closer than 1 1 m to 
the core during the reconstruction. This requirement slightly de- 
creases the efficiency of the array close to tanks and at the same 
time favors core positions further away from the stations. The large 
peak in the center of the array is probably due to this effect com- 
bined with the fact that the nearest station in that direction is par- 
ticularly close. 


7. Determination of energy spectra 

Using the reconstruction methods and quality cuts described in 
Sections 4 and 6, the shower size spectra shown in Fig. 10 were ob- 
tained. In this analysis, the data were split into three zenith angle 
ranges roughly equidistant in seed, defined as: 

Qi = [0°,30°], Q 2 — [30°,40°], n 3 = [40", 46“]. (8) 

A steepening of the spectral slope is visible at log(Si 2 5 /VEM) = 0.5 
and a possible flattening at about log(Si 25 /VEM) = 1.4. To deter- 
mine the energy spectrum from measured data, these S 12 5 spectra 
were unfolded. Unfolding was performed for each zenith angle 
range independently. 

7.3. General method 

For the unfolding procedure the response of the detector to a 
primary particle of mass M, energy £ 0 , zenith angle 0, azimuth <j>, 
and core position ( x c ,y c ) has to be determined from simulation. 
In this analysis we consider only an unfolding of energies. Within 
each zenith angle range, we average over the dependencies on ze- 
nith, azimuth, and core position. The response of the detector is the 
probability of measuring a shower size Si 25 given a primary energy 
E 0 and mass M in a certain zenith range Gl k . 



Fig. 10. Left: Reconstructed shower size spectra in three zenith angle bins. The 
energy spectrum was derived from these spectra using an unfolding method as 
described in Section 7. On the right, the same spectra are shown, weighted with 
Sjj. In this representation the knee can be seen at logS 12 5 ps 0.5 and there are 
indications of a flattening above logSi 2 5 « 1.4. 
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In a discrete formulation we define a response matrix R which 
relates the bin contents N)(i = 1, m) of a measured Si 25 spec- 

trum with the bin contents NJQ' = 1, . . . ,n) of a primary energy 
spectrum for a fixed zenith range Q k : 

N f = R^N]. (9) 

The response matrix elements R® are defined as acceptance 
integrals 


Em dEofdQfdA ± d> M (Eo)p^(S i U5 IE 0 ) 
ZMf K dE of ak dnf Acm dA ± ® M (E 0 ) 


The model flux O m (£o) of nuclei with mass M weighted by their 
acceptance function 

Pm( s 'i25I £ o) =p(Si25,fik|Eo,x c ,y c ,e,4>;!Vf) (11) 

is integrated over primary energy bin E/ 0 , the angles 0 and <j>, and 
area A projected on a plane perpendicular to the particle direction. 
It is summed over all mass components M that contribute to the as- 
sumed composition model. Rp is normalized to the flux integrated 
over bin j in E 0 , solid angle Cl k , and fiducial area A cut . The function p M 
is the probability of an event with mass M and kinematical variables 
(E 0 .x c ,y c , 0, 4>) to be reconstructed with shower size Sj 25 in bin i and 
zenith angle 0 in the range Q k , and to pass all cuts listed in Section 6. 
Thus, Rp for a given primary energy bin j, is the ratio between num- 
ber of events measured in Si 25 bin i and zenith bin k, that pass all 
cuts, and the true number of events in that energy bin j and zenith 
bin k inside the fiducial area. Since the E 0 bins of Rp are indepen- 
dent, the total flux model only affects weighting of events within 
one bin, but not neighboring bins. The flux normalization in R p can- 
cels out, and the dependence on the spectral index of the flux model 
is small (see also Section 8.10). The integrals in Eq. (10) were deter- 
mined numerically using the Monte Carlo method. 

With the normalisation to the full flux integral the response 
matrix has the following normalisation properties (we drop the 
superscript k for zenith range): 

E R 'j = £ j’ E R « = 1 - 02) 

< i 


That means, for a given energy bin j the sum of the probabilities to 
be detected in any signal bin is the efficiency for a given Si 25 bin i 
the probability to belong to any energy E 0 is unity. The efficiency 
depends on the energies, the core position (x c ,y c ) and the angles. 
However, in the analysis we will integrate over core positions and 
azimuth angle. The determination of Ej is described in Eq. (14). 

To obtain the primary energy spectrum from the measured sig- 
nals the matrix Eq. (9) would have to be inverted: 

Nj = (* ')/*• ( 13 ) 

In order to avoid explicit inversion of the badly conditioned matrix 
R, an iterative unfolding method according to Ref. [26] was used. 


7.2. Evaluation of response matrices, efficiencies and resolutions 

Fig. 1 1 shows the response matrix for simulated proton and iron 
primaries in the interval Cii of smallest zenith angles based on an 
E 3 spectrum. In each bin the colour code represents the probabil- 
ity that an event with energy E 0 yields a signal Si 25 . The binning 
uses a logarithmic scale. 

For computational purposes and to smooth fluctuations in the 
simulated response matrix, the logSi 25 projections of each log£ 0 
bin j were fitted by a normal distribution function yielding the 
mean value (logS 125 )j and standard deviation <Ti ogS j. The y 2 proba- 
bility distribution of these fits is almost flat, with a peak at 0 



w 




on 



Fig. 11. Response matrix: shower size S 125 distribution as a function of primary 
energy for proton (left) and iron (right) showers with zenith angles up to 30° 
simulated using the hadronic interaction model SIBYLL. The crosses give the mean 
value and spread (RMS) of the distribution in each energy bin. The matrices have 
been smoothed as described in Section 7.2 and Appendix A. 


caused by cut-off distributions in the threshold region (see also 
discussion of the energy resolution below). The normalisation £j 
was calculated as the ratio between the sum of Monte Carlo event 
weights in the final sample and the sum of weights of events gen- 
erated inside the fiducial area defined in Section 6: 


% = 


EilTw,- 


' J Eil7 w . 


(14) 


Due to migration of shower cores from outside the fiducial area, this 
quantity can become larger than unity. The energy dependence of 
the parameters (logSi 25 ), cr logS , and s was then fitted by empirical 
functions (see A). These functions were used to smooth statistical 
fluctuations in the response matrix and to extrapolate the range 
of simulations to higher energies in order to avoid potential artifacts 
that might be introduced by cutting the spectra off at 100 PeV. 

The mean values and standard deviations are indicated in 
Fig. 1 1 by the points with vertical bars. The average shower sizes 
(logSi 25 (E 0 )) of proton showers for the three different zenith angu- 
lar intervals Q k are shown in Fig. 12(a). Since the shower maximum 
lies above the detector throughout the covered energy range, 
showers from larger zenith angles are more strongly attenuated 
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Fig. 12. (a) Mean shower size as a function of energy for proton showers of various 
zenith angles, (b) Shower size ratio between proton and iron showers. The ratio 
increases for larger zenith angles because the attenuation of iron showers is 
stronger than for proton showers. Both figures are based on an E 3 flux. 


by the atmosphere and thus have a smaller shower size. In 
Fig. 12(b), these points are compared to the mean values for iron. 
The response matrices for proton and iron are very different: on 
average, iron showers have their first interaction at larger height 
leading to a larger shower age than for protons. Iron showers yield 
a smaller average signal than proton showers with the same pri- 
mary energy. The difference between proton and iron increases 
at larger zenith angles. This zenith angle dependence has been 
exploited to test the consistency of our data with models for the 
mass composition, as will be discussed in Section 9. It also means 
that the systematic error due to primary composition will be larger 
for more inclined showers. Therefore, we restricted the final results 
to events with 0 < 30°. 

Fig. 13 shows the efficiencies e obtained in the log£ 0 bins, which 
are the normalisations of the normal distributions of logS^s 
belonging to this bin, for protons and iron nuclei comparing all ze- 
nith angle intervals. The lines are fits to Eq. (A.3). Due to the 
requirement that the station with the largest signal is not on the 
border of the array (see Section 6), peak efficiencies were signifi- 
cantly below 100%. The maximum efficiencies in the three zenith 
angle ranges Cl k correspond to the following effective areas: 

£2, :/W = (1.051 ± 0.013) -10 5 m 2 
n 2 : A ef f = (0.900 ± 0.019) • 10 5 m 2 
n 3 : A eff = (0.803 ± 0.012) • 10 s m 2 . 



Fig. 13. Total efficiency for proton showers in different zenith angle ranges as a 
function of energy. The solid lines are fits of function (A.3), and the dashed lines are 
the results of the corresponding fits for pure iron. 


The decrease of the effective area is a geometrical effect as the pro- 
jection of the fiducial area on the shower plane scales with cos 0. 
Within statistical uncertainties the same values were obtained for 
iron primaries. 

The spread of reconstructed energies (see Fig. 14) has been 
determined by transforming the logSi 2 5 distribution for a given 
E 0 back onto the IogE 0 axis. Above the threshold (between 1 and 
3 PeV depending on zenith angle and primary mass assumption) 
this spread is a measure of the energy resolution. It improves with 
increasing energy, reaching values between 0.04 and 0.12 in 
log(£ 0 ) at 100 PeV, corresponding to a resolution a e /E between 
9% and 23%. The drop below the energy threshold is an artifact of 
the procedure: In the threshold region the trigger condition biases 
toward upward fluctuating showers. This is one reason why this 
region is excluded from the final spectrum. This resolution only 
covers the statistical fluctuations, systematic uncertainties are dis- 
cussed later in Section 8. 

The response matrices obtained by this method depend on the 
primary composition assumption, as well as the hadronic interac- 
tion models and the parametrization of the South Pole atmosphere 
assumed in the simulation. Response matrices were generated for 



Fig. 14. Spread of reconstructed energies as a function of primary energy 
determined as described in the text for proton showers in different zenith angle 
ranges. The lines are fits according to Eq. (A.2) in order to guide the eye. 
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pure proton, pure iron, and for the two-component model de- 
scribed in Section 5.3. 


7.3. Unfolding 


Eq. (9) was solved using an iterative unfolding method based on 
Bayes’ theorem described in Ref. [26], which takes into account the 
total efficiency e and migration due to the fluctuations <7i og (S)- Sim- 
ply inverting the response matrix R would lead to unnatural fluc- 
tuations in the result. 

Starting from a prior distribution P (k) (E 0 j ) (j = 1, . . . ,n) in the 
kth iteration, the inverse of the response matrix R 1 is constructed 
by inverting P(Si 2 5,i|£oj) = Ry using Bayes’ theorem: 


P (fc) (£ 0j |S 12 5,i) 


P(S 125 ,ilEoj)P lk> (E 0 J) 
E/(Si25,#o/)P (k) (£o,.) ' 


(15) 


Then, an estimate of the energy spectrum N e S k \ j = 1 n, is ob- 
tained from the shower size spectrum N s jt i = 1 n: 


NJ 


e(k) 


= lj2NiP {k \E 0 j\Sn5,i)- 


(16) 


In the last step of the iteration, P {k) (E 0 j) is replaced by 

" 7 ) 

As initial prior, P (0) (£ 0 j) ~ Eo 3 was chosen. 

After each iteration, the unfolded spectrum was folded with the 
response matrix, Nf k> = £jRyNj (k \ and compared to the measured 
shower size spectrum. A convergence criterion was then defined 
using the change in y 2 between Nf k) and the measured shower size 
spectrum N 5 between two iterations k and k + 1, as in [27]: 


Ay 2 {k,k+\) = y 2 {N s{k) ,N s ) - y 2 (N s{k+ P,N s ). (18) 


This quantity decreases monotonically during the iteration process. 
However, at A y 2 (k. k + 1) = 0 the unfolding would be equivalent to 
simply inverting R[ k> and the unfolded spectrum would fluctuate 
unnaturally. To avoid this, the iteration was terminated once 
A y 2 (k, k- i-l) fell below a certain value A y 2 erm . The value of this limit 
was determined beforehand using a simple toy simulation in which 
a known spectrum was folded with the response matrix and then, 
after adding statistical fluctuations, unfolded again. In every itera- 
tion step of this unfolding procedure, the unfolded spectrum was 
compared to the known true spectrum. Finally, A y 2 elm = 1.1 was 
chosen where the agreement with the true spectrum was best on 
average. 

The error bars on the unfolded spectrum were determined by 
varying the shower size spectra within their statistical errors and 
repeating the unfolding. This was repeated n = 3000 times and 
the statistical errors in bin j were determined by comparing each 
unfolding result N e|,I) to the average result (N e ): 


(of ) 2 


1 

n- 1 



(19) 


Similarly, bin-to-bin correlations were obtained: 

cov(i,j) = - (J n) (ivf) - (AO,). (20) 

k= 1 


It was verified with a simple toy model that this algorithm correctly 
reproduces a true input spectrum, which was folded with the detec- 
tor response and that the error determination is correct [10]. 


7.4. Correction for snow 

Snow can accumulate at any time on top of the IceTop tanks, 
but a manual measurement of the snow height is only possible 
during the austral summer and therefore is done only once every 
year. The detector simulation took the snow depths measured in 
January, 2007, into account. Data, on the other hand, were taken 
between June and October, 2007, when more snow had accumu- 
lated. In order to estimate the effect of this difference, the detector 
response to proton showers with primary energies of 1 PeV, 10 PeV 
and 30 PeV and zenith angles 0°, 30° and 40° was simulated assum- 
ing once the snow heights measured in January 2007 and once 
those measured in January 2008. In January 2007 the average snow 
depth on top of IceTop tanks was 20.5 cm, while in January 2008 
the average height on top of the same tanks was 53.2 cm. Assum- 
ing constant increase in snow depth and proportionality between 
A log S125 and snow depth, shower sizes in August, 2007, were esti- 
mated. This leads to the following zenith angle dependent energy 
corrections relative to the simulations based on the January 2007 
snow height measurement, which were applied to all unfolded en- 
ergy spectra. Within the statistical uncertainties, no energy depen- 
dence could be observed: 

fi, : Alog(£/PeV) = 0.0368 ± 0.0009, 

Q .2 '. A log(£/PeV) =0.0440 ±0.0013, (21) 

Q 3 : Alog(£/PeV) = 0.0513 ± 0.0008. 


8. Systematic uncertainties 

All systematic errors are summarized in Table 2. In the follow- 
ing, details about the determination of the uncertainties of the en- 
ergy determination and the flux measurement will be given. 

8.1. Snow height 

To estimate the systematic error due to the energy correction 
for snow described in Section 7.4, snow accumulation was as- 
sumed proportional to wind speed. The numbers obtained in this 
way were compared to those assuming constant growth of the 
snow depth (see above). The result of this comparison was used 
as an estimate of the systematic error on energy determination 
due to snow height. 

8.2. Variations of the atmosphere 

As discussed in Section 2.4, variations of the atmosphere affect 
the observed shower sizes. The influence of two parameters of the 
atmosphere has been studied in a data driven way: the total over- 
burden X 0 , and the altitude profile dX,,(/i)/dh. 

First, the days of data taking were ordered according to the total 
atmospheric overburden X 0 . Then the 50 days with the highest and 
the 50 days with the lowest overburden were selected from the to- 
tal of 1 53 days. The average overburdens during these periods 
were X| 0W = 679 g/cm 2 and X high = 700 g/cm 2 , yielding a differ- 
ence of AX = 21 g/cm 2 . From the data taken during these days 
shower size spectra were created for each zenith range O k . 

By comparing the shower size spectra obtained in the two peri- 
ods, the dependence of S125 with atmospheric overburden was de- 
rived. The RMS variation of the total atmospheric overburden 
between June 1 and October 31, 2007, of a Xv = 9-86 g/cm 2 was 
used to estimate the systematic error on the energy determination 
due to atmospheric overburden variations. With the given statisti- 
cal precision, an energy dependence of this variation could not be 
observed. 
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Table 2 

Summary of systematic uncertainties of the energy and flux determination in the three zenith angle intervals fi, , fi 2 . and fi 3 . The individual points are explained in the text. All 
numbers are given for the full energy range, unless otherwise stated. The main source of systematic uncertainty, the unknown primaiy composition, will be treated separately in 
the next section. 


Uncertainty 

o° < e < 30° 


30° < 9 < 40° 


40° ^ 9 < 46° 


Energy (%) 

Flux (%) 

Energy (%) 

Flux (%) 

Energy (%) 

Flux (%) 

Snow height 

0.4 


0.4 


0.4 


Overburden variation 

0.26 


1.9 


3.0 


Atmosphere profile variation 

2.5 


1.8 


1.1 


Atmosphere model 

0.9 


1.1 


0.6 


MC Calibration 

1.5 


1.5 


1.5 


PMT saturation, £ 0 < 10 PeV 

0.5 


0.5 


0.5 


PMT saturation, £ 0 > 1 0 PeV 

<2.5 


<2.5 


<2.5 


Droop 

1.5 


1.5 


1.5 


Calibration stability 

3.0 


3.0 


3.0 


Interaction model, £o < 10 PeV 

2.1 


4.3 


2.0 


Interaction model, £ 0 > 10 PeV 

3.3 


5.5 


4.1 


Shower attenuation 

3.0 


3.0 


3.0 


Flux model 

0.7 


1.0 


1.0 


(logS) and <r logS i 25 

0.7 


1.2 


0.8 


Cut passing rates 


3.0 


3.0 


3.0 

Efficiency 


0.9 


1.6 


1.2 

Unfolding procedure 


1.6 


3.4 


5.2 

Total: E 0 s: lOPeV (%) 

6.0 

3.5 

7.2 

4.8 

6.3 

6.1 

E 0 > lOPeV (%) 

6.9 

3.5 

8.3 

4.8 

7.6 

6.1 


In contrast to total overburden the altitude profile of the atmo- 
sphere at South Pole undergoes a clear annual cycle. To study the 
effect of varying the atmospheric profile on air shower measure- 
ments the data taking period was divided into a period of very 
dense atmosphere (July 25th to October 10th) and one when the 
atmosphere was less dense (remaining days between June 1st 
and July 24th and between October 11th and October 31st). 
Shower size spectra were extracted from the data taken in these 
two periods and by comparing those spectra, an additional system- 
atic error due to the atmospheric profile variation was derived. 

8.3. Atmosphere model in simulation 

The CORSIKA simulations used a model of the South Pole atmo- 
sphere. A systematic uncertainty arises from the choice of model 
since it does not exactly match the average atmosphere during 
the data taking period. To estimate this error on the energy scale 
simulations two different atmosphere parametrizations were com- 
pared. CORSIKA atmosphere model 12 (July 1, 1997), which was 
used in the unfolding procedure, has a total overburden of 
692.9 g/cm 2 and atmosphere model 13 (October 1, 1997) has a to- 
tal overburden of 704.4 g/cm 2 . Averaging the difference in logSi 2 5 
for proton showers between the two simulations above E 0 = 1 PeV 
the systematic error due to the difference of the simulated over- 
burden and the average overburden in data was determined. 

8.4. Calibration 

Systematic uncertainties due to calibration can arise for two 
reasons: variations of the calibration constants between calibra- 
tion runs, and a discrepancy between the calibration of the exper- 
iment and the detector simulation. 

The first point was addressed by studying the variation of the 
number of photoelectrons corresponding to 1 VEM between cali- 
bration runs. From this variation the systematic uncertainty on 
the energy reconstruction due to the tank calibration was esti- 
mated to be 3.0%. 

The simulated tanks were calibrated using the same procedure 
as for the real tanks, as described in Section 5.2. The conversion 
factor between Cherenkov photons and photo electrons resulting 
from this calibration has a statistical uncertainty of 1.5%, which 
was included as a systematic error on the energy. 


8.5. Droop 

The toroid used to decouple the PMT from the signal capture 
electronics introduces a significant droop effect (see Section 2.3), 
which was not corrected for in the analysis. Not correcting for 
droop is not a source of systematic uncertainty in itself if it is done 
consistently in data and simulation. However, discrepancies in the 
way the droop effect is simulated in the detector Monte Carlo, may 
lead to undesired systematic effects. In order to quantify these ef- 
fects, the effect of a droop correction algorithm on the recorded 
charges was compared between data and simulation. From this 
comparison a systematic error on the energy determination of 
1.5% was derived. 

8.6. PMT saturation 

Inaccuracies in the simulated saturation behaviour of the PMT 
could introduce systematic uncertainties on the energy determina- 
tion mostly at high energies. In simulation, saturation sets in at 
higher charges than in the experiment. In order to estimate the ef- 
fect of this discrepancy on the energy spectrum, an artificial, 
charge-based saturation function was applied to the simulated 
charges to bring the simulated charge spectrum into agreement 
with experimental data. Then, the simulated showers were repro- 
cessed, and the change in log(Si 25 /VEM) was used to estimate the 
systematic error on the energy. For primary energies below 10 PeV, 
the systematic error due to the difference in saturation behaviour 
is less than 0.5%. Above 10 PeV it increases exponentially to a value 
of 2.5% at 100 PeV. 

8.7. Cut efficiencies 

Differences in the effects of quality cuts described in Section 6 
when applied to experimental and simulated data lead to a sys- 
tematic uncertainty on the efficiency and consequently on the flux 
normalization. Passing rates of all cuts for data and Monte Carlo 
events above threshold are listed in Table 1. There is a relative dif- 
ference of 3.0% between data and two-component Monte Carlo in 
the total cumulative passing rate, which is included in the system- 
atic uncertainty on the flux. No significant zenith dependence on 
the difference between data and simulation was found. 
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Fig. 15. Shower size ratio for the two-component assumption between SIBYLL and 
QGSJET and EPOS based Monte Carlo simulations for the two component primary 
assumption and showers with zenith angles up to 30°. 

8.8. Interaction model 

Small simulation datasets of proton and iron showers created 
using the high energy hadronic interaction models QGSJET-II and 
EPOS 1.99 in addition to SIBYLL were used to estimate the system- 
atic uncertainty due to the modeling of hadronic interactions. 
Fig. 1 5 shows the shower size ratio between SIBYLL and the alter- 
native simulations as a function of primary energy for the two- 
component primary composition assumption and zenith angles 
up to 30°. Simulations with SIBYLL seem to yield systematically 
smaller shower sizes, and the same observation was made for more 
inclined showers. 

The systematic error derived in this way is purely based on a 
comparison of the three interaction models. All of these models 
have different known strengths and weaknesses in their descrip- 
tion of the underlying physics. Additionally, they all include 
extrapolations of cross-sections and multiplicity distributions to 
energy ranges not accessible by current collider experiments 
which are relevant in the first few cosmic ray interactions. Thus, 
there is an unknown systematic error in case the range of hadron- 
ization models does not cover the true behavior. 


8.9. Shower attenuation 

Since this analysis relies on shower attenuation to extract infor- 
mation about the primary mass composition, we want to ensure 
that the observed attenuation of the measured shower sizes is 
compatible with simulation lying between proton and iron attenu- 
ation. The method of constant intensity cuts [28,29] was used to 
determine shower attenuation. Spectra of N(> S 125 ) in three zenith 
angle ranges were created by integrating the spectra in Fig. 10 left 
for both experimental and simulated data. Assuming an isotropic 
flux of cosmic rays, shower sizes S 125 at different zenith angles that 
result in an equal number of events N(> S 125 ) correspond to the 
same primary energy for any constant mass composition. Shower 
attenuation, i.e. the zenith dependence of S ]2 5 for a given primary 
energy, was measured by determining the values of S 125 in the 
three zenith angle ranges that correspond to six different values 
of N(> S 125 ). These six constant intensity cuts correspond to six dif- 
ferent primary energies. In simulation, intensity levels were chosen 
such that the shower size of proton in the most vertical zenith 
range was the same as in experimental data and for iron cuts were 
made at the same primary energy as for protons. This way no 
assumption on the absolute normalization of the spectrum nor 
an absolute energy scale in experimental data are needed. Fig. 16 
shows the difference between sizes of inclined showers and those 
in the most vertical zenith range. Attenuation in experimental data 
is between the values of proton and iron or consistent with proton. 
As discussed in Section 7.2, iron showers are more strongly atten- 
uated than proton showers. Thus, considering the primary mass 
dependence of shower attenuation the description of attenuation 
in the simulation is consistent with experimental data at the 
reconstruction level. 

However, the predicted shower attenuation differs between ha- 
dronic interaction models. This leads to a zenith dependent sys- 
tematic uncertainty on the change of observed flux for a given 
primary energy with mass. Mass dependence of shower attenua- 
tion shown in Fig. 12(b) for the SIBYLL hadronic interaction model 
was compared to corresponding results obtained with the hadronic 
interaction model EPOS. From the observed zenith dependences of 
the differences between both models a systematic error of the 
shower attenuation of 3% was estimated. Since data from only 
three zenith ranges is available a zenith dependence of this error 
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Fig. 16. Shower attenuation obtained from six constant intensity cuts applied to experimental and simulated shower size spectra as described in Section 8.9. The 
corresponding primary energies in simulation are given in each panel. 


54 


R. Abbasi et al. / Astroparticle Physics 44 (20 13) 40-58 


could not be observed. Note, that this uncertainty may already be 
partially included in the discussion of the systematic error absolute 
energy scale in the previous subsection. 

5.10. Response matrix 

Limited Monte Carlo statistics introduce uncertainties into the 
response matrix. Assuming the efficiency is constant above the 
threshold, the flux error induced by uncertainties of the detector 
response can be estimated by the fit error on c 0 in Eq. (A.3). The 
uncertainties on the parameters a 0 and b 0 in Eqs. (A.1) and (A.2) 
translate to an uncertainty on the energy in the unfolding process. 
These statistical uncertainties on the response matrix were also in- 
cluded in the systematic error of the final result. 

Additionally, the flux model used in the simulation also 
influences the response matrix. A harder spectrum leads to larger 
average shower sizes in an energy bin than a softer one. Simula- 
tions based on an E~ 2 flux and an fT 4 flux were compared with 
the standard simulation which assumes a power law of fT 3 . Above 
the threshold the resulting difference in shower size appears to be 
independent of primary energy. The differences in shower size be- 
tween the two extreme spectral indices were used as an estimate 
of the systematic error on energy scale due to the assumed flux 
model. The resulting systematic uncertainty on the energy spec- 
trum is less than 1%. 

8. 11. Unfolding procedure 

Two parameters besides the response matrix influence the re- 
sult of the unfolding: the termination criterion A^ ax and the prior 
distribution P 0 . Varying the termination criterion, lead to a varia- 
tion of the total flux, which was included as a systematic error. 

In addition, varying the spectral index of the initial prior P 0 be- 
tween -2.5 and -3.5, a variation of the total flux of about 2% was 
observed. Below the knee region around 3 to 4 PeV, the spectral in- 
dex seems to depend on the prior (in the most inclined zenith 
interval even up to 10 PeV. Varying the prior lead to a variation 
of the spectral index below the knee in the most vertical zenith 
band by ±0.01, and in the most inclined zenith range by ±0.025. 
At higher energies variations appear to be purely statistical. 

8.12. Summary of systematic errors 

Systematic uncertainties are summarized in Table 2. The total 
systematic uncertainty was determined by quadratically adding 
the individual contributions. The error on the determination of 
the primary energy in the most vertical zenith angle range is 
6.0% below E 0 = 10 PeV, and 6.9% above. Main contributions are 
the calibration stability (3.0%), atmosphere (2.7% in total), and 
the hadronic interaction model (2.1 %/3.3%). Furthermore, a flux 
uncertainty of 3.5% is caused by differences in cut efficiencies be- 
tween data and Monte Carlo, the efficiency calculation in Monte 
Carlo, and the termination criterion and seed in the unfolding pro- 
cedure. However, the main source of systematic error is the un- 
known primary composition, which will be discussed in detail in 
the next section. 

9. Energy spectrum 

Fig. 17 shows energy spectra for three zenith angular intervals 
unfolded under three assumptions on the mass composition: all- 
proton, all-iron and the two-component model [24] explained in 
Section 5.3. The lower end of the energy range of each spectrum 
was selected where the efficiency according to Eq. (A.3) reached 
90% of the maximum value. The threshold was determined 



log(E/PeV) 



Fig. 17. Resulting flux measured with IceTop, weighted with E 27 . The reconstruc- 
tion was done using three different composition assumptions as described in the 
text: (a) pure proton, (b) pure iron, and (c) Glasstetter's two-component model. In 
each case, the data were divided into three different zenith angle bands equidistant 
in sec(0). Based on the assumption of an isotropic flux, the three individual spectra 
should agree. The boxes indicate the systematic errors. 

individually for each zenith interval and primary composition 
assumption. That way the threshold region is excluded and the 
efficiency can be assumed almost constant. Based on the energy 
resolution, a binning of 10 bins per decade was chosen. 

The systematic error bands were constructed from the numbers 
in Table 2 by defining error boxes around each data point as fol- 
lows: The quadratic sum of flux error and the energy error 
weighted with E 2 7 was used as a vertical error bar. The horizontal 
error bar is given by the error on energy determination. The error 


R. Abbasi et al./Astroparticle Physics 44 (2013) 40-58 


55 


bands were then constructed by connecting the edges of these 
boxes. 

In Fig. 12(b) it was shown that the difference in shower size be- 
tween simulated proton and iron showers increases with zenith 
angle due to the increasing slant depth in the atmosphere, which 
has a different effect for the different masses: iron showers are 
attenuated more strongly with increasing slant depth than proton 
showers. Since the cosmic-ray flux is isotropic to a few per thou- 
sand [30,31] the flux measured in different zenith angular intervals 
has to be the same. 

In case of the pure proton assumption (Fig. 17(a)) a good agree- 
ment between the three spectra is observed. Assuming pure iron 
(Fig. 17(b)), the individual spectra for the three different zenith 
bands clearly disagree at low energies while they start to converge 
towards higher energies. Agreement of the three spectra in case of 
the two-component model (Fig. 17(c)) is good at low and high 
energies. In the intermediate energy range there is some deviation 
between the spectrum obtained from steepest zenith angle range 
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Fig. 18. All particle spectra obtained with IceTop from air showers with zenith 
angles up to 30° under three different composition assumptions: pure proton, the 
two-component model, and a mixture of 30% proton and 70% iron. 


and the other two spectra. However, they are still consistent when 
considering systematic uncertainties. 

Using a y 2 comparison of fluxes in each bin of the spectra from 
the three zenith angle ranges, pure iron could be excluded at a 
>99% confidence level below 24 PeV. This comparison took into ac- 
count both statistical and systematic errors. The latter were treated 
in a conservative way by assuming no correlations between them 
for the different zenith angle intervals. Using the same comparison 
and various mixtures of proton and iron, up to 70% of iron cannot 
be excluded at any energy. A mixture of protons with 70% iron cor- 
respond to a mean logarithmic mass of ln(A) « 2.8. These results 
are in agreement with previous measurements of the composition 
of cosmic rays, which have established a light composition around 
1 PeV and an increasing mass above the knee [32]. 

In Fig. 18, the results obtained in the steepest zenith angle range 
Dt with three primary composition assumptions are compared: 
pure proton, the two-component model, and 70% iron. Only the 
most vertical zenith angle range was chosen, because the differ- 
ence in size for showers initiated by different primaries is smallest 
in this zenith interval, as seen in Fig. 12(b), and because systematic 
uncertainties are smallest in this range. Because the difference in 
shower size between proton and iron decreases toward higher 
energies, the spectrum obtained under the 70% iron assumption 
is softer than the proton-based result. While the composition mod- 
el has a sizable influence on the measured all-particle flux below 
10 PeV, the difference between the two extreme assumptions of 
pure proton and 70% iron almost disappears above 30 PeV. 

As the final result the cosmic ray spectrum is given for the two- 
component model, which represents a realistic mixture of proton 
and iron and is in good agreement with data. Fig. 19 shows the re- 
sult, without the systematic error bands, in comparison to a selec- 
tion of other experiments [33-47], Within systematic errors, our 
results are in good agreement with these previous measurements. 
An analysis of coincident events in IceTop and IceCube [48] using a 
different dataset acquired with a different experimental configura- 
tion resulted in a somewhat lower total flux. However, the results 
are compatible within systematic errors. 

Table 3 lists the measured fluxes including statistical and 
systematic errors. The systematic errors on the flux have been cal- 
culated by transforming the systematic error on energy into a flux 
error based on the local spectral index y : A F/F = yA E/E. This was 
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Fig. 19. The all-particle cosmic ray energy spectrum obtained from the analysis of IceTop data of events with zenith angles up to 30° compared to a selection of other 
experimental results [33-47], 
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Table 3 

All-particle cosmic ray energy spectra measured by the IceTop air shower array for 
the two-component primary composition assumptions using the hadronic interaction 
model SIBYLL2.1. Systematic errors from Table 2 and due to composition are listed 
separately. 


Energy (10 6 GeV) 

dF/dE ± stat ± syst ± comp (GeV 1 m -2 s -1 sr -1 ) 
Two-component assumption 

1.94 

(5.612 ± 0.018 ± 1.2/o|) x 1(T 13 

2.44 

(2.974 ±0.012 ± 0 . 6 / 03 ) x 10 “ 13 

3.07 

(1.589 ± 0.008 ± 0.4/q J 7 ) x 10 13 

3.86 

(8.30 ±0.06 ±1.8/] ;]) x 10 ~ 14 

4.86 

(4.22 ±0.04 ± 1 . 0 / 05 ) x VV 14 

6.12 

(2.098 ± 0.023 ± O. 5 / 023 ) x 10~ 14 

7.71 

(1.021 ±0.015 ± 0 . 24 / 0 ( 12 ) x 10~ 14 

9.70 

(4.92 ± 0.10 ± 1.2 ± 0.5) x 10 - 15 

12.21 

(2.38 ± 0.06 ± 0.6/q J 5 ) x 10 -’ 5 

15.38 

(1.18±0.04±0.29/of ) x 10 - 15 

19.36 

(5.58 ± 0.23 ± 1 ,4/o g 1 ) x 10 ~ 16 

24.37 

(2.66 ± 0.15 ± O. 6 / 021 ) x 10 ~ 16 

30.68 

(1.51 ±0.10 ± O. 4 / 007 5 ) x 10~ 16 

38.62 

(7.5 ± 0.7 ± 1 . 9 / 05 °) x 10 - 17 

48.62 

(3.72 ± 0.4 ± 0.9/g?J°) x 10 17 

61.21 

(1.97 ±0.25 ±0.5/oog°) x 10 -’ 7 

77.06 

(1.05 ±0.17 ±0.27/ooi) x 10 ~ 17 

97.01 

(4.7 ± 1.0 ± 1.4 ± 0.1) x 10 ~ 18 


added quadratically to the systematic error on the flux. The error 
due to the unknown composition is given by the difference be- 
tween the two-component result and the 70% iron and pure proton 
spectra, respectively. 

The spectrum has been fitted with function (6) transformed 
such that the intensity I knee at energy E knee is used as a reference. 
In the fit, statistical errors and bin-to-bin correlations according 
to Eqs. (19) and (20) were used. The results are listed in Table 4. 

The systematic uncertainty of the knee energy £ knee is the sys- 
tematic error on energy determination at that primary energy as 
given in Section 8. The systematic error of J knee has been obtained 
by quadratically adding the systematic error on the flux determi- 
nation and the systematic energy error transformed into a flux 
uncertainty based on the local spectral index. In order to determine 
systematic errors on y, and y 2 , the fit was repeated using the sys- 
tematic errors of the data points as statistical errors. The errors due 
to primary composition were determined by the range of parame- 
ter values obtained when fitting the pure proton and the 70% iron 
results in the same way. The sharpness parameter e of the knee 
was included to obtain an unbiased fit, but is not very well con- 
strained by the data. Its value indicates a relatively sharp knee, 
and in fact we cannot distiguish between a smooth and an infi- 
nitely sharp knee (e — > oo). Fixing s to very large values resulted 
in a j 2 that is not significantly worse. 

Above about 22 PeV a possible flattening of the spectrum can be 
observed independent of primary composition assumption. This 

Table 4 

Fit parameters of the cosmic-ray energy spectrum according to function ( 6 ). 
Systematic errors and uncertainties due to primary composition ("comp”) were 
derived as described in the text. 

Parameter Best fit 

fknee/10 “ 7 m - 2 s - 1 sr - 1 2.38 ±0.23(stat) ±0.5(syst) + 1 - 42 _ 0 . 97 (comp) 

Fknee/PeV 4.32 ±0.22(stat) ±0.26(syst) +° 3S _, ^(comp) 

y 1 -2.759 ±0.015(stat) ±0.5(syst) + 026 _ 0011 (comp) 

y 2 -3.107 ±0.016(stat) ±0.4(syst) + 0 03 _ 008 (comp) 

s 9 ±3(stat) 

X 2/N d{ 19.4/13 


Table 5 

Parameters of the flattening of the spectrum at high 
energy. The errors given are only statistical. 

Ebreak/PeV 23 ±5 

y 3 -2.87 ±0.09 

r7N df 7 . 1/11 


feature is also visible in the measured shower size spectra (see 
Fig. 10). The position of this flattening does not significantly de- 
pend on the primary composition assumption because proton 
and iron showers with 0 < 30° have almost the same shower size 
in this energy range (see Fig. 12(b)). In order to test its statistical 
significance, the spectra were fitted with function (6) plus an addi- 
tional hard break at £ brMk with spectral index y 3 . The goodness of fit 
improves to % 2 /N d f = 7.1/11, which corresponds to a significance 
of 3.2 standard deviations. This, however, does not include system- 
atic errors. The parameters of the flattening are listed in Table 5. 

While a pure power law above the knee cannot be excluded, a 
flattening of the spectrum is preferred by the data. Such a behavior 
has recently been reported by KASCADE-Grande [49]. In addition, 
they observed a steepening of the spectrum just below 10 17 eV, 
which could not be seen in this analysis due to the limited energy 
range. Several models of the cosmic ray energy spectrum predict 
features above the knee similar to the one observed [2]. For in- 
stance, assuming a pure rigidity dependence of the knee there 
can be “gaps” or "concavities” between the knees of two primaries 
if intermediate nuclei are insufficiently abundant to fill in this 
range. Alternatively, a second component of galactic cosmic rays 
could also create the observed behavior. 

10. Summary 

We have derived the all-particle cosmic ray energy spectrum in 
the energy range between 1 PeV and 100 PeV from data taken be- 
tween June and October 2007 with the 26-station configuration of 
the IceTop air shower array at South Pole. 

Using the air shower simulation package CORSIKA with the 
high-energy hadronic interaction model SIBYLL2.1 the relation be- 
tween shower size S125 and primary energy, as well as the detec- 
tion efficiency and energy resolution were determined. Three 
different assumption on the primary mass composition were used 
as input: pure proton, pure iron and a simple two-component 
model [24]. Based on these results, the shower size spectra ob- 
tained in three zenith angle ranges were unfolded with a Bayesian 
unfolding algorithm to obtain energy spectra. 

In case of pure proton and the two-component model, it was 
found that the spectra obtained in the different zenith angle ranges 
were in good agreement. In the pure iron case, on the other hand, a 
strong disagreement between the three spectra was observed at low 
energies. Since one can safely assume that cosmic rays are isotropic 
in the given energy range, the spectra in all three zenith angle ranges 
should be the same. With this assumption, we concluded that pure 
iron primaries can be excluded below energies of 24 PeV. 

We showed that the attenuation of air showers with increasing 
zenith angle bears exploitable information about the chemical 
composition of cosmic rays. Nevertheless, the main source of sys- 
tematic error in the reconstruction of the energy spectrum still re- 
mains the primary mass composition. The systematic error due to 
the choice of a hadronic interaction model is relatively small in this 
analysis because most air shower signals are dominated by the 
electromagnetic component of an air shower, which is relatively 
well understood. For the final result, only the spectra obtained 
from the most vertical zenith angle range, 0° ^ 6 < 30°, were con- 
sidered because in this range the dependence on composition and 
systematic errors are smallest. 


R. Abbasi et al./Astroparticle Physics 44 (2013) 40-58 


57 


By fitting the energy spectrum with function (6), the knee posi- 
tion was determined at 4.3 PeV with spectral indices of -2.76 be- 
low and -3.1 1 above. Around an energy of 22 PeV an indication of 
a flattening of the cosmic ray spectrum to an index of about -2.85 
was observed at the 3u level. 

Since the completion of IceTop and IceCube in 201 1 , the array is 
three times larger than the configuration used in this analysis. 
With this larger array, statistics and containment of high-energy 
showers will be much better, allowing to extend the analysis to 
higher energies. The main strength of IceTop, however, is the pos- 
sibility to measure air showers at the surface in coincidence with 
high-energy muons penetrating deep enough into the ice to trigger 
IceCube. The ratio between the two measurements is highly sensi- 
tive to the mass of the primary particle. 
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Appendix A. Parametrization of the response matrix 


In order to mitigate the effects of statistical fluctuations in the 
unfolding procedure, the response matrices described in Section 7.2 
were separated into mean logarithmic shower size (logS^s), reso- 
lution <7io g s, and efficiency e. There dependences on x = logE p were 
then fitted by empirical functions: 

(logS 125 )(x) = a 0 +x 

ln / exp(ciix) + exp (a 2 + a 3 x + a 4 x 2 ) \ 

V l+exp(a 2 ) )' 


<ho g s(x) 


bo(l + exp (b 3 b 4 )) + exp(-bi)(exp(-b 2 x) - 1) 
1 + exp (— £>3 (x - b 4 )) 


and 


e(x) 


Co 

1 +exp (-c 1 (x-c 2 )+c 3 (x-c 4 )' ! ) 

Co 

l+exp(-c 1 (x-c 2 » 
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X ^ c 4 
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